JSDM models test
Data description
Data was simulated using VirtualCommunity code.
Simulated data contains 20 data sets.
In this file we have compared the 3 models HMSC,GJAM, JSDM on 21 simulated Dataset.\ We have compared the ability to recover the true interactions used for modelling this datasets, by comparing estimated correlation matrix and matrix representing true interactions.
There is a separate file, with the functions needed to fit each model and study the convergence and more parameters are presented.
## Environment filtering 5 species
JSDM
data<-sim_data$EnvEvenSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
Y_cor<-cor(data$Y)
to_prec<-function(m){
n<-dim(m)[1]
Tau_n<-matrix(nrow=n, ncol=n)
for (j in 1:n) {
for (k in 1:n){
Tau_n[j, k] <- -m[j, k]/sqrt((m[j,j]*m[k,k]))
}
}
return(Tau_n)
}
#Tau_n<-matrix(nrow=dim(model$mean$Tau)[1], ncol=dim(model$mean$Tau)[1])
Tau_n<-to_prec(me5$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
cols = colorRampPalette(c("dark blue","white","red"))
col2 <- colorRampPalette(c("#4393C3", "#2166AC", "#053061",
"#FDDBC7", "#FFFFFF", "#D1E5F0", "#92C5DE",
"#67001F", "#B2182B", "#D6604D", "#F4A582"))
corrplot(Y_cor, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(me5$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(me5$mean$EnvRho*(!me5$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(me5$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(me5$mean$Rho*(!me5$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(Tau_n, diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau")
corrplot(Tau_n*(!me5$overlap0$Tau), diag = FALSE, order ="original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")
j_metric_e5<-metrics_jsdm(me5)
#cat("Success rate ")
#j_metric_e5$success_env
cat(sprintf("Success rate: %s\n", metrics_jsdm(me5)$success_env))## Success rate: 0.9mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =as.matrix(me5$mean$Rho*(!me5$overlap0$Rho)))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =as.matrix(me5$mean$Rho*(!me5$overlap0$Rho))) 
GJAM
############################################################Functions
makeSymm <- function(m) {
m[upper.tri(m)] <- t(m)[upper.tri(m)]
return(m)
}
expandSigma_rmd <- function(sigma, S){
ss <- diag(S)
ss[lower.tri(ss,diag=T)] <- sigma
ss[upper.tri(ss)] <- t(ss)[upper.tri(ss)]
ss
}
convert_to_m<-function(ar){
d <-floor((sqrt(length(ar)*8+1)-1)/2)
C <- matrix(0,d,d)
i.lwr <- which(lower.tri(C, diag = TRUE), arr.ind=TRUE)
C[i.lwr] <- ar
C<-makeSymm(C)
return(t(C))
}
metrics_gjam<-function(Rho, comp=NULL,fac=NULL,only_env=T){
if(!is.null(comp) & !is.null(fac)){
fac_comp <- fac-comp
TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac_comp[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac_comp[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac_comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac_comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
TP_comp=TP_comp+1
if(fac_comp[i,j]<0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
wrong=wrong+1
}
}
}else{
if(!is.null(comp)){
TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
TP_fac <- FN_fac<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(comp[i,j]>0 & Rho[i,j] <0)
TP_comp=TP_comp+1
if(comp[i,j]>0 & Rho[i,j] >0)
wrong=wrong+1
if(comp[i,j]>0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if(comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
}
if(!is.null(fac)){
TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
TP_comp <- FN_comp<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac[i,j]>0 & Rho[i,j] <0)
wrong = wrong+1
if(fac[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
}
}
if(only_env){
FP <- TN <- 0
TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):nrow(Rho)){
if(Rho[i,j] != 0)
FP=FP+1
if(Rho[i,j] == 0)
TN=TN+1
}
}
}
success_env<-TN/(TN+FP)
if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
"wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
}
fit_gjam<-function(data, it=2500,burn=500 , name="./gjam_models/temp.rda",interact=diag(ncol(data$Y))){
#setup parameters
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
xdata<-as.data.frame(data$X[,-1])
colnames(xdata)<- c("env","env2")
ydata<-as.data.frame(data$Y)
###formula
formula<-as.formula( ~env+ env2)
ml <- list(ng = it, burnin = burn, typeNames = 'PA')
####fit
mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
save(mod_gjam1, file = name)
#summary(mod_gjam1)
hist(effectiveSize(mod_gjam1$chains$sgibbs), main="ess(beta)",lwd=2,col=gray(.6))
hist(effectiveSize(mod_gjam1$chains$bgibbs), main="ess(beta)",lwd=2,col=gray(.6))
# postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
# postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
pH<-convert_to_m(postH)
pL<-convert_to_m(postL)
R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
sgibbs<-mod_gjam1$chains$sgibbs
tau<-array(NA,dim=c(data$J,data$J,1000))
for(j in 1:1000){
ss <- expandSigma_rmd(sgibbs[j,], S = data$J)
si <- solve(ss)
tau[,,j] <- -cov2cor(si)
}
# sgibbs<-mod_gjam1$chains$sgibbs
# tau<-array(NA,dim=c(data$J,data$J,it-burn))
# for(j in (burn:it)){
# ss <- expandSigma_rmd(sgibbs[j,], S = 5)
# si <- solve(ss)
# tau[,,j-burn] <- -cov2cor(si)
# }
tau_mean<-apply(tau,c(1,2), mean)
tau_HI<-apply(tau,c(1,2),quantile,0.95)
tau_LO<-apply(tau,c(1,2),quantile,0.05)
Tau_sign<-tau_mean*(!(tau_HI>0 & tau_LO<0))
# par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
# corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("Correlation cor(Y)")
# corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("R")
# corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("E matrix")
# corrplot(Tau_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("Tau signif")
# corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("R signif")
# corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("True interactions")
R_sign
}
###this function only loads the model and return the R and T significant
load_gjam<-function(data,it=2500,burn=500,name="./gjam_models/temp.rda",interact=diag(ncol(data$Y)),autoC=F){
#setup parameters
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
gj_mod<-load_object(name)
#summary(gj_mod)
S<-ncol(data$Y)
ind <-which(array(diag(S)[lower.tri(diag(S),diag=T)])==1)
par(mfrow=c(1,3),oma = c(1, 1, 1, 1))
hist(effectiveSize(gj_mod$chains$sgibbs[,-ind]), main="sigma non diagonal elements",lwd=2,col=gray(.6))
hist(effectiveSize(gj_mod$chains$sgibbs[,ind]), main="sigma diagonal elements",lwd=2,col=gray(.6))
hist(effectiveSize(gj_mod$chains$bgibbs), main="ess(beta)",lwd=2,col=gray(.6))
if(autoC){
plot(acfplot(mcmc(gj_mod$chains$bgibbs)))
plot(acfplot(mcmc(gj_mod$chains$sgibbs)))
}
postH<-apply(gj_mod$chains$sgibbs[burn:dim(gj_mod$chains$sgibbs)[1],], 2, quantile,0.95)
postL<-apply(gj_mod$chains$sgibbs[burn:dim(gj_mod$chains$sgibbs)[1],], 2, quantile,0.05)
pH<-convert_to_m(postH)
pL<-convert_to_m(postL)
R_sign<-cov2cor(gj_mod$parameters$sigMu)*(!(pH>0 & pL<0))
sgibbs<-gj_mod$chains$sgibbs
tau<-array(NA,dim=c(data$J,data$J,it-burn))
for(j in (burn:it)){
ss <- expandSigma_rmd(sgibbs[j,], S = data$J)
si <- solve(ss)
tau[,,j-burn] <- -cov2cor(si)
}
tau_mean<-apply(tau,c(1,2), mean)
tau_HI<-apply(tau,c(1,2),quantile,0.95)
tau_LO<-apply(tau,c(1,2),quantile,0.05)
Tau_sign<-tau_mean*(!(tau_HI>0 & tau_LO<0))
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(data$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(gj_mod$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(gj_mod$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("E matrix")
corrplot(Tau_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Tau signif")
corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R signif")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
return(list(Rho_sign=R_sign,Tau_sign=Tau_sign))
}
########################################################################to setup chains parameters
data<-sim_data$EnvEvenSp5
#Rho_sign<-fit_gjam(data,2500,500,name="./gjam_models/gjam5env.rda")
mod_gjam<-load_gjam(data,name="./gjam_models/gjam5env.rda", autoC=T)
g_metric_e5<-metrics_gjam(mod_gjam$Rho_sign)
cat(" ")cat(sprintf("Success rate: %s\n", metrics_gjam(mod_gjam$Rho_sign)$success_env))## Success rate: 0.5mod_list_Rho<-list.append(mod_list_Rho,gjam=mod_gjam$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=mod_gjam$Tau)HMSC
#setwd("~/Tesi/Code/Ecology-models-master/simcoms-master/ExampleFiles")
fit_hmsc<-function(data,label="Fit",nsamples = 1000,nchains=2,name="./HMmodels/hmtemp.rda" ){
if (label=="Fit"){
Y_data = subset(data, select = -env)
ns<- ncol(Y_data)
np <- nrow(Y_data)
X<-scale(poly(data$env[1:np], 2))
colnames(X)<-c("env","env2")
studyDesign = data.frame(sample = as.factor(1:np))
rL = HmscRandomLevel(units = studyDesign$sample)
m = Hmsc(Y=as.matrix(Y_data), XData=as.data.frame(X), XFormula=~env+env2, distr="probit",
studyDesign = studyDesign, ranLevels = list(sample = rL))
m = sampleMcmc(m, nsamples, thin=10, adaptNf=c(200,200), transient=500,nChains=nchains ,verbose=F)
save(m, file = name)
return(m)
}
if (label=="Load"){
return(load_object(name))
}
}
data<-sim_data$EnvEvenSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5env.rda" )
#hm_mod<-load_object("./HMmodels/hm5env.rda")Convergence:
hm_conv<-function(mod, autoC=F){
codaList = convertToCodaObject(mod)
#convergence histograms
hist(effectiveSize(codaList$Beta), main="ess(beta)",lwd=2,col=gray(.6))
hist(gelman.diag(codaList$Beta,multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(beta)")
if(autoC){
plot(acfplot(codaList$Beta))
plot(acfplot(codaList$Omega[[1]]))
}
hist(effectiveSize(codaList$Omega[[1]]), main="ess(omega)",lwd=2,col=gray(.6))
hist(gelman.diag(codaList$Omega[[1]], multivariate=FALSE)$psrf,lwd=2,col=gray(.6), main="psrf(omega)")
}
hm_conv(mod=hm_mod,autoC=TRUE )
Study of interactions
metrics_hmsc<-function(Rho, comp=NULL,fac=NULL,only_env=T){
#{
if(!is.null(comp) & !is.null(fac)){
fac_comp <- fac-comp
TP_comp <- FN_comp <-TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac_comp[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac_comp[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac_comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac_comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
if(fac_comp[i,j]< 0 & Rho[i,j] < 0)
TP_comp=TP_comp+1
if(fac_comp[i,j]<0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if((fac_comp[i,j]<0 & Rho[i,j] > 0) | (fac_comp[i,j]>0 & Rho[i,j] < 0))
wrong=wrong+1
}
}
}else{
if(!is.null(comp)){
TP_comp <- FN_comp <- FP <- TN <- wrong <- 0
TP_fac <- FN_fac<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(comp[i,j]>0 & Rho[i,j] <0)
TP_comp=TP_comp+1
if(comp[i,j]>0 & Rho[i,j] >0)
wrong=wrong+1
if(comp[i,j]>0 & Rho[i,j] == 0)
FN_comp=FN_comp+1
if(comp[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(comp[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
# }
# }
}
if(!is.null(fac)){
TP_fac <- FN_fac <- FP <- TN <- wrong <- 0
TP_comp <- FN_comp<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):ncol(Rho)){
if(fac[i,j]>0 & Rho[i,j] >0)
TP_fac=TP_fac+1
if(fac[i,j]>0 & Rho[i,j] <0)
wrong = wrong+1
if(fac[i,j]>0 & Rho[i,j] == 0)
FN_fac=FN_fac+1
if(fac[i,j]==0 & Rho[i,j] == 0)
TN=TN+1
if(fac[i,j]==0 & Rho[i,j] != 0)
FP=FP+1
}
}
}
}
if(only_env){
FP <- TN <- 0
TP_fac <- FN_fac<-TP_comp <- FN_comp<-wrong<-NULL
for(i in 1:(nrow(Rho)-1)){
for(j in (i+1):nrow(Rho)){
if(Rho[i,j] != 0)
FP=FP+1
if(Rho[i,j] == 0)
TN=TN+1
}
}
}
success_env<-TN/(TN+FP)
if(!is.null(TP_comp)){success_comp <- TP_comp/(TP_comp+FN_comp)}else{ success_comp = NULL}
if(!is.null(TP_fac)) {success_fac <- TP_fac/(TP_fac+FN_fac)}else{ success_fac = NULL}
list("FP"=FP,"TN"=TN,"TP_comp"=TP_comp,"FN_comp"=FN_comp,"TP_fac"=TP_fac,"FN_fac"=FN_fac,
"wrong"=wrong,"success_env"=success_env,"success_comp"=success_comp,"success_fac"=success_fac)
}
hm_inter<-function(mod, nchains=2,nsamples = 1000, interact=diag(ns)){
getOmega = function(a,r=1)
return(crossprod(a$Lambda[[r]]))
ns<-mod$ns
postOmega1 = array(unlist(lapply(mod$postList[[1]],getOmega)),c(ns,ns,mod$samples))
postOmega2 = array(unlist(lapply(mod$postList[[2]],getOmega)),c(ns,ns,mod$samples))
postOmega<-abind(postOmega1,postOmega2,along=3)
postOmegaMean = apply(postOmega,c(1,2),mean)
postOmegaUp=apply(postOmega,c(1,2),quantile,0.95)
postOmegaLo=apply(postOmega,c(1,2),quantile,0.05)
postR<-array(dim=c(ns,ns,nchains*nsamples))
postT<-array(dim=c(ns,ns,nchains*nsamples))
for(i in 1:dim(postOmega)[3]){
postR[,,i]<-stats::cov2cor(postOmega[,,i])
postT[,,i]<- -stats::cov2cor(solve(postOmega[,,i]+diag(ns)))
}
postRMean = apply(postR,c(1,2),mean)
postRUp=apply(postR,c(1,2),quantile,0.95)
postRLo=apply(postR,c(1,2),quantile,0.05)
postTMean = apply(postT,c(1,2),mean)
postTUp=apply(postT,c(1,2),quantile,0.95)
postTLo=apply(postT,c(1,2),quantile,0.05)
Toplot_R<-postRMean*(!(postRUp>0 & postRLo<0))
Toplot_T<-postTMean*(!(postTUp>0 & postTLo<0))
# Omegacor<- computeAssociations(m)
# supportLevel<- 0.95
# toPlot<- ((Omegacor[[1]]$support>supportLevel)+ (Omegacor[[1]]$support<(1-supportLevel))>0)*Omegacor[[1]]$mean
# corrplot(toPlot, method="color", col=colorRampPalette(c("blue", "white", "red"))(200))
par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
corrplot(cor(hm_mod$Y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(postRMean, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R")
corrplot(Toplot_R, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Plot only non zero value")
corrplot(Toplot_T, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Partial correlation matrix")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
return(list(Rho_sign=Toplot_R,Tau_sign=Toplot_T))
}
hm_mod_R<-hm_inter(hm_mod,nchains=2,nsamples = 3000)
h_metric_e5<-metrics_hmsc(hm_mod_R$Rho_sign)
# cat("Success rate ")
# h_metric_e5$success_env
cat(" ")
cat(sprintf("Success rate: %s\n", metrics_hmsc(hm_mod_R$Rho_sign)$success_env))
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R JSDM ")
corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R GJAM")
corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("R HMSC")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
}
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(5))
ALL4<-function(jsdm_mod,gjam_mod,hmsc_mod,interact=diag(5)){
par(mfrow=c(2,2),oma = c(1, 1, 1, 1))
#corrplot(jsdm_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
#title("Partial correlation JSDM ")
corrplot(gjam_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Partial correlation GJAM")
corrplot(hmsc_mod, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Partial correlation HMSC")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
}
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(5))
Environmental filtering 10 species
JSDM
me10 <- load_object("model-2019-04-10-08-26-20.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me10)## Summary for model '/tmp/RtmpKix4lZ/file40e957831178'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 683.734 minutes at time 2019-04-09 21:02:35.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.jsdm_conv(me10)## Maximum Rhat value for Beta: 1.5602
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.7953
me10$mcmc.info[1:7]## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000j_metric_e10<-metrics_jsdm(me10)
cat("Success rate ")## Success ratej_metric_e10$success_env## [1] 0.9333333data<-sim_data$EnvEvenSp10
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
plot_cor_jsdm<-function(mod,y,interact=diag(ncol(y))){
par(mfrow=c(2,4),oma = c(3, 1, 2, 1))
corrplot(cor(y), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Correlation cor(Y)")
corrplot(mod$mean$EnvRho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho")
corrplot(mod$mean$EnvRho*(!mod$overlap0$EnvRho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("EnvRho signif")
corrplot(mod$mean$Rho, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho")
corrplot(mod$mean$Rho*(!mod$overlap0$Rho), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("Rho signif")
corrplot(interact, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
title("True interactions")
}
#plot_cor_jsdm(me10,data$Y)
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =me10$mean$Rho*(!me10$overlap0$Rho))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =me10$mean$Rho*(!me10$overlap0$Rho)) Gjam
data<-sim_data$EnvEvenSp10
#Rho_sign<-fit_gjam(data,3000,500,"./gjam_models/gjam10env.rda")
gjam_mod<-load_gjam(data,3000,500,name="./gjam_models/gjam10env.rda")
#gje10<-load_object("./gjam_models/gjam10env.rda")
g_metric_e10<-metrics_gjam(gjam_mod$Rho_sign)
cat("Success rate ")## Success rateg_metric_e10$success_env## [1] 0.3111111#to check posterior density of s in Sigma
#gje10<-load_object("./gjam_models/gjam10env.rda")
#plot(density(gje10$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$EnvEvenSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10env.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod)
h_metric_e10<-metrics_hmsc(hm_mod_R$Rho_sign)
# cat("Success rate ")
# h_metric_e10$success_env
# h_metric_e10$success_fac
cat(" ")cat(sprintf("Success rate: %s\n", h_metric_e10$success_env))## Success rate: 1cat(sprintf("Success rate: %s\n", h_metric_e10$success_fac))
mod_list_Rho<-list.append(mod_list_Rho, hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, hmsc=hm_mod_R$Tau_sign) 
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(10))
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(10))
Environmental filtering 20 species
JSDM
me20 <- load_object("model-2019-04-11-19-06-02.rda")
#load("model-2019-04-09-19-02-16.rda")
summary(me20)## Summary for model '/tmp/RtmpKix4lZ/file40e966e51ba7'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2079.661 minutes at time 2019-04-10 08:26:21.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.me20$mcmc.info[1:7]## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000jsdm_conv(me20)## Maximum Rhat value for Beta: 1.2174
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3236
j_metric_e20<-metrics_jsdm(me20)
cat("Success rate ")## Success ratej_metric_e20$success_env## [1] 0.9894737data<-sim_data$EnvEvenSp20
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
plot_cor_jsdm(me20,data$Y)
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =me20$mean$Rho*(!me20$overlap0$Rho))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =me20$mean$Rho*(!me20$overlap0$Rho)) 
Gjam
data<-sim_data$EnvEvenSp20
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
gjam_mod<-load_gjam(data,5000,500,name="./gjam_models/gjam20env.rda")
#gje20<-load_object("./gjam_models/gjam20env.rda")
g_metric_e20<-metrics_gjam(gjam_mod$Rho_sign)
cat("Success rate ")## Success rate#g_metric_e20$success_env
cat(" ")cat(sprintf("Success rate env: %s\n", round(metrics_gjam(gjam_mod$Rho_sign)$success_env,3)))## Success rate env: 0.358#to check posterior density of s in Sigma
#gje20<-load_object("./gjam_models/gjam20env.rda")
#plot(density(gje20$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) Gjam dimension reduction
#
# data<-sim_data$EnvEvenSp20
#
# #fit_gjam(data,5000,500,"./gjam_models/gjam20env.rda")
#
# #load_gjam(data,name="./gjam_models/gjam20env.rda")
# #gje20<-load_object("./gjam_models/gjam20env.rda")
#
#
# #to check posterior density of s in Sigma
# #gje20<-load_object("./gjam_models/gjam20env.rda")
# #plot(density(gje20$chains$sgibbs[,4]))
# data <- list(
# Y = subset(data, select = -env),
# X = cbind(1, scale(poly(data$env, 2))),
# covx = cov(cbind(1, scale(poly(data$env, 2)))),
# K = 3,
# J = ncol(data) - 1,
# n = nrow(data),
# I = diag(ncol(data) - 1),
# df = ncol(data)
# )
# xdata<-as.data.frame(data$X[,-1])
# colnames(xdata)<- c("env","env2")
# ydata<-as.data.frame(data$Y)
# ###formula
# rl <- list(r = 8, N = 20)
# formula<-as.formula( ~env+ env2)
# ml <- list(ng = 2500, burnin = 500, typeNames = 'PA', reductList = rl)
# ####fit
#
#
# mod_gjam1 <- gjam(formula, xdata = xdata, ydata = ydata, modelList = ml)
# save(mod_gjam1, file = "./gjam_models/gjam20env_dr.rda")
# summary(mod_gjam1)
#
# Tau <- solve(mod_gjam1$parameters$sigMu)
# Tau_n = to_prec(Tau)
#
# #postH<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.95)
# #postL<-apply(mod_gjam1$chains$sgibbs, 2, quantile,0.05)
#
#
# #pH<-convert_to_m(postH)
# #pL<-convert_to_m(postL)
#
# #R_sign<-cov2cor(mod_gjam1$parameters$sigMu)*(!(pH>0 & pL<0))
#
# g_metric_e20_dim_red<-metrics_gjam(Rho_sign)
# cat("Success rate ")
# g_metric_e20_dim_red$success_env
#
# par(mfrow=c(2,3),oma = c(1, 1, 1, 1))
# corrplot(cor(ydata), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("Correlation cor(Y)")
# corrplot(mod_gjam1$parameters$corMu, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("R")
# corrplot(mod_gjam1$parameters$ematrix, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("E matrix")
# # corrplot(Tau_n, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# # title("Tau")
# # corrplot(R_sign, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", type = "lower")
# # title("R signif")
# corrplot(diag(20), diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color",col=cols(200), type = "lower")
# title("True interactions")
# data<-sim_data$EnvEvenSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20env.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod)
h_metric_e20<-metrics_hmsc(hm_mod_R$Rho_sign)
cat("Success score ")## Success scoreh_metric_e20$success_env## [1] 1cat(" ")cat(sprintf("Success rate non-interacting: %s\n", round(metrics_hmsc(hm_mod_R$Rho_sign)$success_env,3)))## Success rate non-interacting: 1mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,diag(20))
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,diag(20))
Environmental filtering + Facilitation dense 5 species
JSDM
mfd5 <- load_object("model-2019-04-11-19-35-11.rda")
summary(mfd5)## Summary for model '/tmp/RtmpKix4lZ/file40e94e482135'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.13 minutes at time 2019-04-11 19:06:03.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.#mfd5$Rhat
jsdm_conv(mfd5)## Maximum Rhat value for Beta: 1.3081
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1421
mfd5$mcmc.info[1:7]## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000j_metric_facDense5<-metrics_jsdm(mfd5,fac = fac_inter[[4]],only_env = F)
cat("Success rate ")## Success ratej_metric_facDense5$success_env## [1] 1j_metric_facDense5$success_fac## [1] 0data<-sim_data$FacDenseSp5
data <- list(
Y = subset(data, select = -env),
X = cbind(1, scale(poly(data$env, 2))),
covx = cov(cbind(1, scale(poly(data$env, 2)))),
K = 3,
J = ncol(data) - 1,
n = nrow(data),
I = diag(ncol(data) - 1),
df = ncol(data)
)
##########################################################################################
#Tau<-solve
#Tau_n<-to_prec(me10$mean$Tau)
#Tau_k<-Tau_n*(!(model$q97.5$Tau>0 & model$q2.5$Tau<0))
#plot_cor_jsdm(mfd5,data$Y,fac_inter[[4]])
mod_list_Rho<-list()
mod_list_Rho<-list(jsdm =mfd5$mean$Rho*(!mfd5$overlap0$Rho))
mod_list_Tau<-list()
mod_list_Tau<-list(jsdm =mfd5$mean$Rho*(!mfd5$overlap0$Rho)) Gjam
data<-sim_data$FacDenseSp5
#Rho_sign<-fit_gjam(data,10000,2000,"./gjam_models/gjam5f.rda",interact=fac_inter[[4]])
gjam_mod<-load_gjam(data,10000,2000,name="./gjam_models/gjam5f.rda", interact=fac_inter[[4]])
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
g_metric_facDense5<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[4]], only_env = F)
#cat("Success rate ")
#g_metric_facDense5$success_env
#g_metric_facDense5$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5f.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense5$success_env,3)))## Success rate for non-interacting: 0.5cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense5$success_fac,3)))## Success rate for facilitation: 0.5mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign)HMSC
data<-sim_data$FacDenseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fd.rda" )
hm_conv(hm_mod)
hm_inter(hm_mod, interact = fac_inter[[4]])## $Rho_sign
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 0 0
## [3,] 0 0 1 0 0
## [4,] 0 0 0 1 0
## [5,] 0 0 0 0 1
##
## $Tau_sign
## [,1] [,2] [,3] [,4] [,5]
## [1,] -1 0 0 0 0
## [2,] 0 -1 0 0 0
## [3,] 0 0 -1 0 0
## [4,] 0 0 0 -1 0
## [5,] 0 0 0 0 -1hm_mod_R<-hm_inter(hm_mod)
h_metric_facDense5<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[4]],only_env = F)
#cat("Success rate ")
# h_metric_facDense5$success_env
# h_metric_facDense5$success_fac
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense5$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense5$success_fac,3)))## Success rate for facilitation: 0mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[4]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[4]])
Environmental filtering + Facilitation dense 10 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e96d09e31e'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 684.487 minutes at time 2019-04-11 19:35:12.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.2866
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4808
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate## [1] 0.925## [1] 0data<-sim_data$FacDenseSp10
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10fd.rda",interact=fac_inter[[5]])
gjam_mod<- load_gjam(data,5000,500,name="./gjam_models/gjam10fd.rda", interact=fac_inter[[5]])
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
g_metric_facDense10<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[5]], only_env = F)
# cat("Success rate ")
# g_metric_facDense10$success_env
# g_metric_facDense10$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense10$success_env,3)))## Success rate for non-interacting: 0.625cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense10$success_fac,3)))## Success rate for facilitation: 0.8mod_list_Rho<-list.append(mod_list_Rho,gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,gjam=gjam_mod$Tau_sign)HMSC
## Success rate for non-interacting: 1## Success rate for facilitation: 0
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[5]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[5]])
Environmental filtering + Facilitation dense 20 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e91ec9ff9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2062.467 minutes at time 2019-04-12 06:59:42.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.6132
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.508
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate## [1] 0.9833333## [1] 0Gjam
data<-sim_data$FacDenseSp20
#Rho_sign<-fit_gjam(data,10000,1500,"./gjam_models/gjam20fd.rda",interact=fac_inter[[6]])
gjam_mod<-load_gjam(data,10000,1500,name="./gjam_models/gjam20fd.rda", interact=fac_inter[[6]])
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
g_metric_facDense20<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[6]], only_env = F)
# cat("Success rate ")
# g_metric_facDense20$success_env
# g_metric_facDense20$success_fac
#
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facDense20$success_env,3)))## Success rate for non-interacting: 0.406cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facDense20$success_fac,3)))## Success rate for facilitation: 0.5mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacDenseSp20
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm20fd.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[6]])
h_metric_facDense20<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[6]],only_env = F)
# cat("Success rate ")
# h_metric_facDense20$success_env
# h_metric_facDense20$success_fac
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facDense20$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facDense20$success_fac,3)))## Success rate for facilitation: 0mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[6]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[6]])
Environmental filtering + Facilitation sparse 5 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e917d13dae'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.05 minutes at time 2019-04-13 17:22:13.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.373
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4048
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.889## Success rate for facilitation: 1Gjam
data<-sim_data$FacSparseSp5
#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam5fs.rda",interact=fac_inter[[7]])
gjam_mod<-load_gjam(data,2500,500,name="./gjam_models/gjam5fs.rda", interact=fac_inter[[7]])
#gjfd5<-load_object("./gjam_models/gjam5fs.rda")
g_metric_facSparse5<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[7]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse5$success_env
# g_metric_facSparse5$success_fac
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse5$success_env,3)))## Success rate for non-interacting: 0.556cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse5$success_fac,3)))## Success rate for facilitation: 0#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacSparseSp5
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm5fs.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[7]])
h_metric_facSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse5$success_env
# h_metric_facSparse5$success_fac
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse5$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse5$success_fac,3)))## Success rate for facilitation: 0mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[7]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[7]])
Environmental filtering + Facilitation sparse 10 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9113ddc9'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 686.687 minutes at time 2019-04-13 17:51:16.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.5479
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4746
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.93## Success rate for facilitation: 0data<-sim_data$FacSparseSp10
#Rho_sign<-fit_gjam(data,2500,500,"./gjam_models/gjam10fs.rda",interact=fac_inter[[8]])
gjam_mod<-load_gjam(data,name="./gjam_models/gjam10fs.rda", interact=fac_inter[[8]])
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
g_metric_facSparse10<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[8]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse10$success_env
# g_metric_facSparse10$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10fs.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse10$success_env,3)))## Success rate for non-interacting: 0.326cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse10$success_fac,3)))## Success rate for facilitation: 0mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacSparseSp10
hm_mod<-fit_hmsc(data,"Load",name="./HMmodels/hm10fs.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, interact = fac_inter[[8]])
h_metric_facSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[8]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse10$success_env
# h_metric_facSparse10$success_fac
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse10$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse10$success_fac,3)))## Success rate for facilitation: 0mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[8]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[8]])
Environmental filtering + Facilitation sparse 20 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9629c97ad'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2056.316 minutes at time 2019-04-14 05:17:59.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.5653
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.4826
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.978## Success rate for facilitation: 0Gjam
data<-sim_data$FacSparseSp20
#fit_gjam(data,10000,1500,"./gjam_models/gjam20fs.rda",interact=fac_inter[[9]])
gjam_mod <-load_gjam(data,10000,1500,name="./gjam_models/gjam20fs.rda", interact=fac_inter[[9]])
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
g_metric_facSparse20<-metrics_gjam(gjam_mod$Rho_sign,fac=fac_inter[[9]], only_env = F)
# cat("Success rate ")
# g_metric_facSparse20$success_env
# g_metric_facSparse20$success_fac
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_facSparse20$success_env,3)))## Success rate for non-interacting: 0.446cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_facSparse20$success_fac,3)))## Success rate for facilitation: 0.25mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
HMSC
data<-sim_data$FacSparseSp20
hm_mod<-fit_hmsc(data,nsamples=3000, nchains=2,"Load",name="./HMmodels/hm20fs.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = fac_inter[[9]])
h_metric_facSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,fac=fac_inter[[9]],only_env = F)
# cat("Success rate ")
# h_metric_facSparse20$success_env
# h_metric_facSparse20$success_fac
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_facSparse20$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_facSparse20$success_fac,3)))## Success rate for facilitation: 0mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,fac_inter[[9]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,fac_inter[[9]])
Environmental filtering + Competition dense 5 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9127bbb96'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.308 minutes at time 2019-04-15 15:34:20.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## Successful convergence based on Rhat values (all < 1.1).## Maximum Rhat value for Beta: 1.0324
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.0313
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate## [1] 1## [1] 1## Success rate for non-interacting: 1## Success rate for facilitation: 0Gjam
data<-sim_data$CompDenseSp5
#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam5cmpd.rda",interact=comp_inter[[10]])
gjam_mod<-load_gjam(data,10000,1000,name="./gjam_models/gjam5cmpd.rda", interact=(-1)*comp_inter[[10]])
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
g_metric_compDense5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[10]], only_env = F)
# cat("Success rate ")
# g_metric_compDense5$success_env
# g_metric_compDense5$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense5$success_env,3)))## Success rate for non-interacting: 0.625cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_compDense5$success_comp,3)))## Success rate for facilitation: 1#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5cmpd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) data<-sim_data$CompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hm5cmpd.rda" )
#hm_mod<-fit_hmsc(data,"Fit",nsamples=3000, nchains=2 )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[10]])
h_metric_compDense5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[10]],only_env = F)
# cat("Success rate ")
# h_metric_compDense5$success_env
# h_metric_compDense5$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense5$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_compDense5$success_comp,3)))## Success rate for facilitation: 1mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[10]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[10]])
Environmental filtering + Competition dense 10 species
## Summary for model '/tmp/RtmpKix4lZ/file40e9559ce9d8'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 695.614 minutes at time 2019-04-15 16:03:39.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.1845
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1771
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.875## Success rate for facilitation: 1Gjam
data<-sim_data$CompDenseSp10
#Rho_sign<-fit_gjam(data,5000,500,"./gjam_models/gjam10cmpd.rda",interact=comp_inter[[11]])
gjam_mod<-load_gjam(data,5000,500,name="./gjam_models/gjam10cmpd.rda", interact= (-1)*comp_inter[[11]])
#gjfd5<-load_object("./gjam_models/gjam10cmpd.rda")
g_metric_compDense10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[11]], only_env = F)
# cat("Success rate ")
# g_metric_compDense10$success_env
# g_metric_compDense10$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense10$success_env,3)))## Success rate for non-interacting: 0.325cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense10$success_comp,3)))## Success rate for competition: 1#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$CompDenseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm10cmpd.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod,nsamples=5000, nchains=2, interact = (-1)*comp_inter[[11]])
h_metric_compDense10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[11]],only_env = F)
# cat("Success rate ")
# h_metric_compDense10$success_env
# h_metric_compDense10$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense10$success_env,3)))## Success rate for non-interacting: 0.875cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense10$success_comp,3)))## Success rate for competition: 1mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[11]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[11]])
Environmental filtering + Competition dense 20 species
## Summary for model '/tmp/RtmpKix4lZ/file40e96843124a'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2061.424 minutes at time 2019-04-16 03:39:17.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.9065
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.6944
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.956## Success rate for competition: 0.7Gjam
data<-sim_data$CompDenseSp20
#Rho_sign<-fit_gjam(data,10000,1000,"./gjam_models/gjam20cmpd.rda",interact= (-1)*comp_inter[[12]])
gjam_mod<-load_gjam(data,10000,1000,name="./gjam_models/gjam20cmpd.rda", interact= (-1)*comp_inter[[12]])
#gjfd5<-load_object("./gjam_models/gjam20cmpd.rda")
g_metric_compDense20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[12]], only_env = F)
# cat("Success rate ")
# g_metric_compDense20$success_env
# g_metric_compDense20$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compDense20$success_env,3)))## Success rate for non-interacting: 0.394cat(sprintf("Success rate for competition: %s\n", round(g_metric_compDense20$success_comp,3)))## Success rate for competition: 0.9#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20fd.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$CompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=5000, nchains=2,name="./HMmodels/hm20cmpd.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=5000, nchains=2,interact = (-1)*comp_inter[[12]])
h_metric_compDense20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[12]],only_env = F)
# cat("Success rate ")
# h_metric_compDense20$success_env
# h_metric_compDense20$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compDense20$success_env,3)))## Success rate for non-interacting: 0.967cat(sprintf("Success rate for competition: %s\n", round(h_metric_compDense20$success_comp,3)))## Success rate for competition: 0.4#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[12]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[12]])
Environmental filtering + Competition sparse 5 species
## Summary for model '/tmp/RtmpKix4lZ/file40e97e225117'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.251 minutes at time 2019-04-17 14:00:46.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.1876
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2038
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 1## Success rate for competition: 1Gjam
data<-sim_data$CompSparseSp5
#Rho_sign<-fit_gjam(data,5000,2000,"./gjam_models/gjam5cmps.rda",interact= (-1)*comp_inter[[13]])
gjam_mod<-load_gjam(data,5000,2000,name="./gjam_models/gjam5cmps.rda", interact= (-1)*comp_inter[[13]])
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
g_metric_compSparse5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[13]], only_env = F)
cat("Success rate ")## Success rateg_metric_compSparse5$success_env## [1] 0.5555556g_metric_compSparse5$success_comp## [1] 1cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse5$success_env,3)))## Success rate for non-interacting: 0.556cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse5$success_comp,3)))## Success rate for competition: 1#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam5cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$CompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm5cmps.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[13]])
h_metric_compSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[13]],only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse5$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse5$success_comp,3)))## Success rate for competition: 1#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[13]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[13]])
Environmental filtering + Competition sparse 10 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e978d641ed'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 685.719 minutes at time 2019-04-17 14:30:01.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.116
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.1753
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate## [1] 0.9534884## [1] 1## Success rate for non-interacting: 0.953## Success rate for competition: 1Gjam
data<-sim_data$CompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam10cmps.rda",interact= (-1)*comp_inter[[14]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjam10cmps.rda", interact= (-1)*comp_inter[[14]])
g_metric_compSparse10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[14]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse10$success_env
# g_metric_compSparse10$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse10$success_env,3)))## Success rate for non-interacting: 0.326cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse10$success_comp,3)))## Success rate for competition: 1#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam10cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) data<-sim_data$CompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm10cmps.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[14]])
h_metric_compSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[14]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse10$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse10$success_comp,3)))
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[14]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[14]])
Environmental filtering + Competition sparse 20 species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9379ee963'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2084.707 minutes at time 2019-04-18 01:55:46.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.2652
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3883
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.968## Success rate for competition: 0Gjam
data<-sim_data$CompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjam20cmps.rda",interact= (-1)*comp_inter[[15]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjam20cmps.rda", interact= (-1)*comp_inter[[15]])
g_metric_compSparse20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[15]], only_env = F)
# cat("Success rate ")
# g_metric_compSparse20$success_env
# g_metric_compSparse20$success_comp
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_compSparse20$success_env,3)))## Success rate for non-interacting: 0.269cat(sprintf("Success rate for competition: %s\n", round(g_metric_compSparse20$success_comp,3)))## Success rate for competition: 0.75#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$CompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hm20cmps.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[15]])
h_metric_compSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[15]],only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_compSparse20$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for competition: %s\n", round(h_metric_compSparse20$success_comp,3)))## Success rate for competition: 0#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc,(-1)*comp_inter[[15]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc,(-1)*comp_inter[[15]])
Environmental filtering + Competition +Facilitation 5 dense species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e93b49ad29'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 28.267 minutes at time 2019-04-19 12:40:31.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.1667
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2086
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.667## Success rate for competition: 1## Success rate for facilitation: 0.5Gjam
data<-sim_data$FacCompDenseSp5
#Rho_sign<-fit_gjam(data,30000,15000,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
gjam_mod<-load_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign3<-fit_gjam(data,5000,1000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign4<-fit_gjam(data,5000,2500,"./gjam_models/gjam5cmp_facd.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign5<-fit_gjam(data,20000,5000,"./gjam_models/gjam5cmp_facd2.rda",interact= (-1)*comp_inter[[16]]+fac_inter[[16]])
#Rho_sign<-load_gjam(data,name="./gjam_models/gjam5cmp_facd.rda", (-1)*comp_inter[[16]]+fac_inter[[16]])
g_metric_FacCompDense5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[16]], fac=fac_inter[[16]], only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense5$success_env,3)))## Success rate for non-interacting: 0.833cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense5$success_comp,3)))## Success rate for competition: 1cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense5$success_fac,3)))## Success rate for facilitation: 1# par(mfrow=c(2,3))
#corrplot((Rho_sign+Rho_sign2)/2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign2, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign3, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign4, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
# corrplot(Rho_sign5, diag = FALSE, order = "original",tl.pos = "ld", tl.cex = 0.5, method = "color", col=cols(200),type = "lower")
#to check posterior density of s in Sigma
#gjfd5<-load_object("./gjam_models/gjam20cmps.rda")
#plot(density(gjfd5$chains$sgibbs[,4]))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacCompDenseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=3000, nchains=2,name="./HMmodels/hmcmp_facd5.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=3000, nchains=2,interact = (-1)*comp_inter[[16]]+ fac_inter[[16]])
h_metric_FacCompDense5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[16]],fac= fac_inter[[16]],only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense5$success_env,3)))## Success rate for non-interacting: 0.833cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense5$success_comp,3)))## Success rate for competition: 0.5cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense5$success_fac,3)))## Success rate for facilitation: 0.5#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[16]]+ fac_inter[[16]])
Environmental filtering + Competition +Facilitation 10 dense species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e93b173410'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 688.628 minutes at time 2019-04-19 13:08:47.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.3458
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3195
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.914## Success rate for competition: 0## Success rate for facilitation: 0Gjam
data<-sim_data$FacCompDenseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd10.rda",interact= (-1)*comp_inter[[17]]+fac_inter[[17]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facd10.rda", interact= (-1)*comp_inter[[17]]+fac_inter[[17]])
g_metric_FacCompDense10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[17]], fac=fac_inter[[17]], only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense10$success_env,3)))## Success rate for non-interacting: 0.286cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense10$success_comp,3)))## Success rate for competition: 0.6cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense10$success_fac,3)))## Success rate for facilitation: 1mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd10.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[17]] + fac_inter[[17]])
h_metric_FacCompDense10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[17]],fac=fac_inter[[17]], only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense10$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense10$success_comp,3)))## Success rate for competition: 0cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense10$success_fac,3)))## Success rate for facilitation: 0#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[17]]+ fac_inter[[17]])
Environmental filtering + Competition +Facilitation 20 dense species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e95eadc1e3'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2073.109 minutes at time 2019-04-20 00:37:26.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.3008
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.2755
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.965## Success rate for competition: 0.333## Success rate for facilitation: 0.333Gjam
data<-sim_data$FacCompDenseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facd20.rda",interact= (-1)*comp_inter[[18]]+fac_inter[[18]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facd20.rda", interact= (-1)*comp_inter[[18]]+fac_inter[[18]])
g_metric_FacCompDense20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[18]], fac=fac_inter[[18]], only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompDense20$success_env,3)))## Success rate for non-interacting: 0.366cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompDense20$success_comp,3)))## Success rate for competition: 0.667cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompDense20$success_fac,3)))## Success rate for facilitation: 0.429mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacCompDenseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facd20.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[18]] + fac_inter[[18]])
h_metric_FacCompDense20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[18]],fac=fac_inter[[18]],only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompDense20$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompDense20$success_comp,3)))## Success rate for competition: 0cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompDense20$success_fac,3)))## Success rate for facilitation: 0#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[18]]+ fac_inter[[18]])
Environmental filtering + Competition +Facilitation 5 sparse species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e9896a5fb'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 29.285 minutes at time 2019-04-21 11:10:35.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.3771
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3486
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.875## Success rate for competition: 1## Success rate for facilitation: 1Gjam
data<-sim_data$FacCompSparseSp5
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs5.rda",interact= (-1)*comp_inter[[19]]+ fac_inter[[19]] )
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs5.rda", interact= (-1)*comp_inter[[19]]+ fac_inter[[19]])
g_metric_FacCompSparse5<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[15]], only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse5$success_env,3)))## Success rate for non-interacting: 0.556cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse5$success_comp,3)))## Success rate for competition: 1cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse5$success_comp,3)))## Success rate for competition: 1mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacCompSparseSp5
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs5.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[19]]+ fac_inter[[19]])
h_metric_FacCompSparse5<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[19]],fac= fac_inter[[19]] ,only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse5$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse5$success_comp,3)))## Success rate for competition: 1cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse5$success_fac,3)))## Success rate for facilitation: 0#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[19]]+ fac_inter[[19]])
Environmental filtering + Competition +Facilitation 10 sparse species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e95f910eba'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 679.081 minutes at time 2019-04-21 11:39:53.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 2.1146
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 3.3048
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.951## Success rate for competition: 0## Success rate for facilitation: 0Gjam
data<-sim_data$FacCompSparseSp10
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs10.rda",interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs10.rda", interact= (-1)*comp_inter[[20]] +fac_inter[[20]])
g_metric_FacCompSparse10<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]], only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse10$success_env,3)))## Success rate for non-interacting: 0.341cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse10$success_comp,3)))## Success rate for competition: 1cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparse10$success_fac,3)))## Success rate for facilitation: 0mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacCompSparseSp10
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs10.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[20]] +fac_inter[[20]])
h_metric_FacCompSparse10<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[20]],fac=fac_inter[[20]],only_env = F)
cat("\n")cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse10$success_env,3)))## Success rate for non-interacting: 1cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse10$success_comp,3)))## Success rate for competition: 0cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse10$success_fac,3)))## Success rate for facilitation: 0#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[20]]+ fac_inter[[20]])
Environmental filtering + Competition +Facilitation 20 sparse species
JSDM
## Summary for model '/tmp/RtmpKix4lZ/file40e978c52358'
## Saved parameters: B Rho EnvRho
## MCMC ran in parallel for 2083.296 minutes at time 2019-04-21 22:58:59.
##
## For each of 5 chains:
## Adaptation: 250000 iterations (possibly insufficient)
## Burn-in: 0 iterations
## Thin rate: 100 iterations
## Total chain length: 270000 iterations
## Posterior sample size: 200 draws
##
## **WARNING** Rhat values indicate convergence failure.## Maximum Rhat value for Beta: 1.2707
## Maximum Rhat value for Rho: NA
## Maximum Rhat value for EnvRho: 1.3875
## $n.chains
## [1] 5
##
## $n.adapt
## [1] 250000 250000 250000 250000 250000
##
## $sufficient.adapt
## [1] FALSE FALSE FALSE FALSE FALSE
##
## $n.iter
## [1] 20000
##
## $n.burnin
## [1] 0
##
## $n.thin
## [1] 100
##
## $n.samples
## [1] 1000## Success rate for non-interacting: 0.984## Success rate for competition: 0## Success rate for facilitation: 0Gjam
data<-sim_data$FacCompSparseSp20
#Rho_sign<-fit_gjam(data,2000,500,"./gjam_models/gjamcmp_facs20.rda",interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
gjam_mod<-load_gjam(data,2000,500,name="./gjam_models/gjamcmp_facs20.rda", interact= (-1)*comp_inter[[21]]+fac_inter[[21]])
g_metric_FacCompSparse20<-metrics_gjam(gjam_mod$Rho_sign,comp=comp_inter[[21]], fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(g_metric_FacCompSparse20$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(g_metric_FacCompSparse20$success_comp,3)))
cat(sprintf("Success rate for facilitation: %s\n", round(g_metric_FacCompSparse20$success_fac,3)))
mod_list_Rho<-list.append(mod_list_Rho, gjam=gjam_mod$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau, gjam=gjam_mod$Tau_sign) HMSC
data<-sim_data$FacCompSparseSp20
hm_mod<-fit_hmsc(data,"Load",nsamples=2000, nchains=2,name="./HMmodels/hmcmp_facs20.rda" )
hm_conv(hm_mod)
hm_mod_R<-hm_inter(hm_mod, nsamples=2000, nchains=2,interact = (-1)*comp_inter[[21]] +fac_inter[[21]])
h_metric_FacCompSparse20<-metrics_hmsc(hm_mod_R$Rho_sign,comp=comp_inter[[21]],fac=fac_inter[[21]],only_env = F)
cat("\n")
cat(sprintf("Success rate for non-interacting: %s\n", round(h_metric_FacCompSparse20$success_env,3)))
cat(sprintf("Success rate for competition: %s\n", round(h_metric_FacCompSparse20$success_comp,3)))
cat(sprintf("Success rate for facilitation: %s\n", round(h_metric_FacCompSparse20$success_fac,3)))
#mod_list<-list.append(mod_list,hmsc=Rho_sign)
mod_list_Rho<-list.append(mod_list_Rho,hmsc=hm_mod_R$Rho_sign)
mod_list_Tau<-list.append(mod_list_Tau,hmsc=hm_mod_R$Tau_sign)
ALL3(mod_list_Rho$jsdm,mod_list_Rho$gjam,mod_list_Rho$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])
ALL4(mod_list_Tau$jsdm,mod_list_Tau$gjam,mod_list_Tau$hmsc, (-1)*comp_inter[[21]]+ fac_inter[[21]])
Graph from paper
Visualization for true interactions
Facilitation
Competition Dense
Competition Dense
Competition + Facilitation Dense
Competition + Facilitation Sparse
